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Abstract 


The steady state simulation of multistage, multicomponent 
separators requires the solution of a set of nonlinear algebraic 
equations. The Newton-Raphson technique is commonly used to solve 
the set of equations. The Napthali-Sandholm method for single 
separation columns gives rise to a block-tridiagonal Jacobian. 
Most of the blocks are highly sparse. A modified Gaussian 
elimination has been developed to solve the jacobian Many 
separation problems occur at high pressures where gamma-phi 
approach will lead to errors. For those cases the Equation of state 
approach has been incorporated .The latest mixing rule given by 
Wong and Sandler has been used with Peng-Robinson equation of 
state to calculate K in this present work. Computer programs have 
been written implementing the modified VLE calculations and 
successfully tested with absorption and distillation problems. 



Chapter : 1 
Introduction 

Separation processes like absorption, extraction, 
distillation etc. are widely used in chemical process 
industries to separate binary or multicomponent mixtures. 
Separation units are often the major factor in plant capital and 
operating costs. Therefore design, modelling and simulation of 
these units are extremely important. In this work, stagewise 
contact separation columns have been considered. 

The modern day computers have drastically reduced the manual 
effort involved in solving separation problems in a rigorous 
manner. Many approaches have been proposed in the literature for 
efficient solution of the nonlinear algebraic equations involving 
mass balance, heat balance and equilibrium relations which arise 
in the modelling of the separation columns. The simultaneous 
convergence method proposed by Naphtali and Sandholm [5] is a 
well established procedure to simulate single separation columns. 

In the Naphtali-Sandholm method, 
the set of nonlinear algebraic equations are solved by the 
Newton-Raphson technique. The iterative solution involves a 
Jacobian, which is a large and sparse matrix. In case of single 
separation columns the Jacobian has a block-tridiagonal 
structure. Most of the constituent blocks of Jacobian are 
highly sparse. The large and sparse structure of the Jacobian can i 
be exploited in solving the linear equations efficiently. The 
Newton-Raphson method may fail to converge to solution if the i 
initial guess is not within the domain of convergence and 
more over, only one solution can be obtained from a particular 



starting point even if there are multiple solutions. These 
limitations of the Newton-Raphson technique can be overcome by the 
use of homotopy continuation methods. These methods convert the 
nonlinear algebraic equations to a set of coupled nonlinear 
ordinary differential equations with initial conditions. These 
equations are solved by predictor - corrector methods which trace 
out the solution path starting from the initial condition. The 
actual solution procedure involve the solution of large and sparse 
matrices for every predictor-corrector step. The structure of 
these matrices are that of the Jacobian of the Newton-Raphson 
method with one additional row and a column at the end. Efficient 
solution of the matrices are a crucial part of the solution 
procedure. 

The initial guess for Newton-Raphson method can be generated 
in various ways. However, there is no single scheme which can 
ensure convergence for all types of separation problems. Recently, 
Venkata Ramana[7] proposed a simple initialisation scheme based 
feed conditions. 

The objective of the present work is : 

To incorporate high pressure phase equilibria based on 
equation of state approach. 

In chapter 2 of the present work, a brief review of the 
literature avaliable on the relevant topics is presented. Chapter 
3 contains the modelling equations and the methods of solution by 
Newton-Raphson procedure using the proposed schemes, for 
single columns . Chapter 4 describes the theoretical aspects of 



Equation of State approach. Finally, results and discussions 
are presented in Chapter 5 . 



Chapter : 2 
Literature Review 


The steady state simulation of a stagewise contact separation 
column involves the determination of the temperatures and 
component flow rates for every stage and at the ends of the column 
under specified conditions. The basic procedure is to solve the 
component mass balances, equilibrium relations and enthalpy 
balances for all the stages. The large number of variables and 
complexity of the equations have led to many different solution 
techniques.’ A good review of the previous work has been presented 
by Venkata Ramana [7]. 

In a classic paper in 1971, Naphtali and Sandholm[5] proposed 
simultaneous convergence scheme for all the stage temperatures and 
component flow rates. They grouped by generating the equations in 
a stage by stage fashion. The Murphree stage efficiencies were 
included in the equlibrium relations. 

The set of nonlinear equations are solved by using the 
Newton-Raphson method. The Jacobian which arises in the Newton 
-Raphson procedure has a large and sparse structure. In case of 
single columns, the Jacobian has a block tri-diagonal structure. 

A brief review of the suggested methods of solving the Jacobian 
is presented below. 

1) Single column : Naphtali and Sandholm used the standard 
Gauss elimination scheme to solve the tridiagonal Jacobian. 
Henley and Seader [3] extended the Thomas algorithm for 



tridiagonal matrices to the case of the block-tridiagonal 
Jacobian. Their procedure, known as the block Thomas algorithm, is 
a popular method to solve the Jacobian for single columns. 

2 . Continuation methods : The Newton-Raphson technique as 

applied to separation problems fail to converge to a solution if 
the initial guess is not in the domain of convergence. Also the 
method produces only one solution from a given starting point. In 
cases where the Newton-Raphson method fails to converge to a 
solution or where multiple solutions are being seeked, the 
continuation methods are employed. Jelinek et al. [9] appear to 
be the first to employ the concept of continuation to separation 
problems. They obtained solutions for multicomponent multistage 
separation with constant flow rate using reflux ratio as the 
continuation parameter. As mentioned by Wayburn and Seader [8], 
attempts have been made to apply continuation to separation 
problems using different homotopy functions, attempting 
stabilisation of the ODE-IVP formed from the algebraic equations 
and also by path parametrisation. They suggested an affine 
homotopy which prevents the solution path doubling back to the 
initial guess value in place of the simple Newton's homotopy. 
Wayburn and Seader [10] obtained multiple solutions for 
interlinked distillation column problems and showed applicability 
of continuation in cases where the Newton-Raphson method failed. A 
recent review of the previous work can be found elsewhere [7]. 
Recently Wong and Sandler [11] proposed a newmixing rule. This new 
mixing rule with Peng-Robinson Equation of State has been used 
to calculate thr equilibrium constants [ K] which is required when 



one uses the Murphree efficiency as the equilibrium relation. 



Chapter : 3 

Modelling of Multistage Multicomponent Separation Columns 
AND Methods of Solution 


In this chapter the modelling equations for multistage 
multicomponent separators and their methods of solution are 
presented. The model proposed by Naphtali and Sandholm [5] for 
single columns has been followed with certain modifications as 
discussed by Venkata Ramana [7]. A method of solution of the 
Jacobian using a modified Gaussian elimination scheme has been 
presented. 

A generalised modified Gaussian elimination scheme has also 
been proposed to handle the different possible configurations of 
the Jacobian. 


3.1 Modelling Equations for Single Columns : 

A counter-current steady state vapour-liquid separator 
consisting of N stages and separating a C component mixture is 
considered. The stages of the separator are numbered from the top 
to bottom. A schematic representation of a general stage j of the 
separator is depicted in Figure-1. As can be seen from Figure -1 
the liquid and vapour leaving the stage are denoted by and 
respectively, the vapour from the lower stage and the liquid from 


the upper stage are denoted by V. and L._ repectively, the stage 

j ■'■i- 3 X 

has a feed f^ and the liquid and vapour side-streams leaving the 
stage are denoted by and respectively. A heat 

input/withdrawal of q^ is also considered. 




F1G.1 A GENERAL STAGE IN A SINGLE COLUMN 



The modelling equations for this general stage are the 


component material balances, equilibrium relations and the 
enthalpy balance. These equations when expressed as discrepancy 
functions appear as follows : 

(a) Component Mass Balance : 

....a...... ^3.1.1) 

where ,l:si^C , & ,lsj:£N 

(b) Equilibrium Relations : 


''j.i 

Q . . = 1 . . 

L. V. 


"^j+l, i 


+ (l-T?. .) 


V. 


D+1 


V 


1/1 


^j,i ^j,i" ^j+l,i 


. (3.1.2) 
.(3.1.3) 


where ,l5iiC , & ,lij£N and t). . is the Murphree 

efficiency for component i at stage j . 

(c) Enthalpy Balance : 


"j = £ E -a+s.) J 


-( 1 +S.) TH. .V. .+rh^ f. . + Q. 

J' j' 1,1 1,1 ^ 

1 ^ j N and 1 i i :£ C 

(3.1.4) 

Sj and Sj are the fractions of the vapour and liquid 
side-streams from the j-th stage respectively. 



The above equations form a set of (2C+1) equations in terms 
of (2C+1) variables for each stage and hence a total of N(2C+1) 
equations in N(2C+l) unknowns for the entire column. The above 
equations can be written in a compact form as 
F(X) = 0 


. (3.1.5) 


where 


F = 

(Fi, F^, . 

‘'j 

(3.1.6) 

F. = 

'”3.1' •• 

,M. p, Q. , .. ,Q. p, E. )'^ 

J Ja-*“ J r ^ J 

(3.1.7) 

X = 

(Xj^, x^. 

'X., ...., X^)^ 

(3.1.8) 

X. = 

3 

'h,i' •• 

'^j,C' ^j,l' * ■ 

(3.1.9) 


As mentioned by Venkata Ramana [7], the equations and 
variables are ordered in such a way so as to get a matrix 
structure whose sparsity can be suitably exploited. 

The Newton-Raphson method can be applied to obtain the 
solution of equation (3.1.5). The correction vector AX is given by 

-111 ^ 


Ax" = 


^SF] 

I 5X J . 


. . . (3.1.10) 


or, AX^ = - ( J "^)^ ...(3.1.11) 

where J is the Jacobian matrix and I is the iteration number. 

For the first iteration, an initial guess value of X, X° has 
. 1 

to be estimated and AX calculated from equation (3.1.11). Then 
the I-th approximation of the variables is computed as 

X^ = X^"^ + e AX^ ... (3.1.12) 

The value of 0 can be chosen in such a way that the 


convergence is either accelerated or dampened. In the usual 
Newton-Raphson technique 0 is taken as one. Henley and Seader 
[3] suggested a Fibonacci search to find an optimum 0 for each 



iteration. Venkata Ramana [7] adopted the Broyden's method to 
calculate 8 at each iteration step and the same has been used in 
the present work. 

3.2 Method of Solution of Jacobian for Single Column : 

The structure of the Jacobian matrix obtained from the 
Naphtali-Sandholm formulation is block-tridiagonal in case of 
single columns. The structure of the Jacobian is given by 



The structures of the A, B and C blocks are given by: 
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where I is a C x C identity matrix, O is a C x C null matrix, O 
c cc c 

a C dimensional null vector, is a C x C nonzero matrix, is 

cc c 

a C dimensional nonzero vector and X is a nonzero number. As can 
be seen, A is a highly sparse matrix, C is also sparse while B has 
few zeros. 

The matrix equation involving the tridiagonal Jacobian 
can be solved by the block Thomas algorithm. However, a Gaussian 
elimination which takes into consideration the highly sparse 
structure of the Jacobian reduces the operation count 
significantly as will be shown later. 

A method of solution for the linearised matrix equation 
J AX = -F has been developed in the present work based on the 
Gaussian elimination technique. The proposed algorithm takes into 
account the sparse structure of the Jacobian and its constituent 
submatrices (blocks A, B and C) to reduce the operation count as 
well as storage requirements as compared to a full matrix or a 
block-tridiagonal matrix with filled up blocks. 

The Gaussian elimination : 

The Gaussian elimination technique is discussed briefly 
now. Consider a general matrix equation M x = b, where M is a 
matrix of order m x m, x and b are m dimensional vectors. Let the 



elements of M, x and b be as follows: 



— 


— 


^11 

^12 

^Im 

M = 

^21 

^22 

^2m 


^ml 

^m2 

• • • • d. 

mm 


X = Lx^ Xj ... x^ ; b = |_b^ bj ... b^Il'^ 

The augmented matrix of the above matrix equation is defined by 
the matrix |~M |b~|. The Gaussian elimination procedure for this 
general matrix equation consist of two parts : 

i) The forward elimination - A series of operations are 
carried out on the equations such that the M matrix is transformed 
into an upper-triangular matrix with diagonal elements as one. 
This is achieved by dividing the i-th row of the augmented matrix 
by the element at location (i,i), i.e., the pivot element 

(normalisation step) . Then the resultant row is multiplied by the 
element at (l,i) of M and subtracted from the 1-th row where i+1 
1 :£ m (reduction step). The above steps are carried out 
successively from i = 1 to i = m, ultimately yielding the 

upper-triangular matrix. 

ii) The back substitution - The last element of the 

augmented matrix after the forward elimination is x . Now, the x.s 

m i 

are solved for with the help of known x^s (i+1 £ 1 ^ m) by moving 
up the rows of the augmented matrix. 



The algorithm consists of the following steps : 


Forward elimination 


1. a^^^ < 



2 . bj^ < 

bi/aii ] 

1:5 i :< m ; i+1 

► 

3. < 

■ ®11 J 

i+1 S k :s m 

4. b,^ < 

‘’k - ^kl ‘>1 


Back substitution 


5- = ’=!» 

in 


6. X. = b. 

■ ? ^ii ""i 
l=i+l 

m-1 2: i 2 : 1 


m 


Modified Gaussian elimination for single column : 

The Gaussian elimination has been suitably modified in 
the present work to exploit the block-tridiagonal structure of 
the Jacobian in case of a single column. Furthermore, the presence 
of zeros and ones in the individual blocks has been exploited to 
reduce the operation count, i.e., number of multiplications and 
divisions performed and also storage requirements. 

Since all the diagonal elements (pivots) of the Jacobian 
are also the diagonal elements of B blocks, the forward 
elimination can be thought of as transforming B blocks to an 
upper-triangular form and reducing all the elements of the blocks 
directly below the B blocks to zeros. This way a block by block 
approach has been formulated. For the j-th B block, the only 
nonzero block directly below it is the block 


because of the 



tridiagonal structure of the Jacobian, and hence we have to 
operate only on the j-th and (j+l)-th rows of blocks for the j-th 
row. Also, there are no nonzero elements directly above and 
hence this block remains unaltered throughout the forward 
elimination of the j-th row and need not be considered. So, the 
blocks which get affected at the j-th step of the elimination are 
Bj , Cj , Fj , ^j+1 * *’j reduces to zero 
at the (j-l)-th step and therefore is not present for forward 
elimination of the j-th row. 

From the above discussion it follows that the operations 
for the j-th row of Jacobian will involve 4C+2 rows (2C+1 rows 
each for the j and (j+l)-th blocks) and 4C+3 columns (2C+1 columns 
corresponding to each of and and one corresponding to the F 
blocks) . The operations can be easily understood by following the 
structure of the Jacobian at different stages during the forward 
elimination as shown in the following discussion. 

Forward substitution : 

For the first B block, the pivots are the elements (1,1) 
to (2C+1,2C+1) . The structure of the blocks which are operated 
upon in the augmented Jacobian is 
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where x represents nonzero numbers. 

For the first C pivots the pivot row is divided by the 
quantity (1+s^) if s^ 0, if s^ = 0 (as shown in the above 
diagram for the sake of simplicity) the divisions are not 
performed. Division is performed only on elements in columns C+i, 
3C+l+i and 4C+3 . As the rest of the elements of the row are zeros, 
they are not operated upon. The elements directly below the pivot 
element are reduced to zeros. Only the rows C+1 to 2C+1, 2C4-i+i 

dnd 4C+2 are to be operated , since all the other rows below the 
pivot are already in the reduced form as can be verified from the 
structure of the individual blocks of the Jacobian. For the 
(2C+l+i)-th row no multiplication is involved as the element in 
the pivot column is unity. Again, only the elements in columns 
C+i, 3C+l+i and 4C+3 are affected for the rows operated upon. 

For pivots (C+1, C+1) to (2C,2C), the pivot row is 
divided by the pivot element, the elements which are divided are 
in columns C+i+1 to 2C+1, 3C+2 to 4C+1 and 4C+3. Other elements of 
the row are zeros and hence not divided. The rows which are 


considered for the reduction step are C+i+1 to 2C+i+i and 4C+2, 
the rest being already in the reduced form because of the 
block structure. Also, only the elements in the columns C+i+1 to 
2C+1, 3C+2 to 4C+1 and 4C+3 are affected. 

For pivot (2C+1,2C+1), the (2C+l)-th row is divided by 
the (2C+1, 2C+1) th element. The columns 2C+2 to 3C+1 are not 
considered, since the elements are zeros. The rows 2C+2 to 3C+1 
and 4C+2 are reduced. Again, rows 3C+2 to 4C+1 are not operated 
upon due to the reduced structure they have in the original 
blocks. Elements corresponding to columns 2C+2 to 3C+1 are left 
out of the operations. 

The above operations complete the first stage and the 
blocks under consideration now have the following structure 
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The elements which have filled in the zeros of the blocks or had 
their values changed are denoted by x’ in the above matrix. 

It is seen that the block becomes upper-triangular 
with the diagonal elements as unity, the block retains its 



structure 

except 

the 

identity 

C X C block which now 

becomes a 

diagonal 

matrix. 

The 

identity 

block remains unchanged 

if s^ is 

equal to 

zero . 

The 

A_ block 

becomes a null matrix 

after the 


operations. The block becomes more filled up, with its second C 
X C diagonal matrix becoming a full matrix and the zero vector 
also becoming filled up. 

A point to be noted here is that the first C columns of 
the block remain unchanged because the first C columns of the 
block are all zeros. Also the (C+l)-th to 2C-th rows of B^ are 
not affected because the corresponding block elements are zeros 
in the original matrix and remain so through out. These two 
features are retained through out the forward elimination. 

The nonzero elements of the blocks B^ and and the 

vector are now stored for the back substitution. 

The blocks B^ , , A^ , B^ and F^ are to be 

considered for the second stage of the forward elimination. 

The steps for pivots (1,1) to (C,C) are similar to those 
for the first row except that the elements in the columns C+l to 
2C+1 have to be operated upon in addition to the columns which are 
operated upon for the first row of blocks. 

For pivots (C+l, C+l) to (2C+1,2C+1), the operations are 
the same as those for row 1. 

The structure of the Jacobian after the second B block 
has been reduced is as shown below 




Only the nonzero elements of the and blocks and the 

vector are stored for back substitution. 

The point to note here is that the structure of the 
block now is same as the structure that the block B^ had at 
the end of the forward elimination of the first row of blocks. The 
implication is that from here onwards the steps of the forward 
elimination will exactly be the same upto block and the 

upper-triangular B blocks, its adjacent C blocks and the 
corresponding F blocks will be stored suitably at the end of the 
operations for each row of blocks. 

For the last B block, i.e. B„ , the blocks involved in 
the operations are the B^^ block itself and the F^^ vector. No C or 
A blocks are to be considered in this case as B^^ does not have any 
adjacent C block or lower A block to it. The operations involved 
are similar to those for the previous blocks, but there are no 
operations for row or column numbers greater than 2C+i except the 
(4C+3)-th column corresponding to F^^ . As before, the 

upper-triangular B block and the F block are stored in this case 



for the back substitution. 


Back substitution : 

The upper-triangular matrix is now easily solved by back 
substitution starting with the last element of the N-th block. For 
any AX(j,i) solved for in the N-th block the values of the solved 
AX(j,k)s are used, k going from i+1 to 2C+1. For AX(j,i), the i-th 
row of the N-th block is considered and AX(j,k)s corresponding to 
the zeros in that row of the upper-triangular matrix are not 
considered for operations. 

For blocks N-1 to 1, the AX values of only the (j+l)-th 
block have to be considered in addition to the j-th block for the 
j-th row of blocks. This is because of the block-tridiagonal 
structure of the Jacobian. For the i-th equation of the j-th 
block, all the (j+l)-th AXs corresponding to the nonzero elements 
of the (j+l)-th block in that row are considered. In addition, the 
AX(j,k)s of the j-th block corresponding to the nonzero elements 
in the i-th row are also to be considered (for i+l:sj<2C+l) . 

Implementation of the algorithm : 

As mentioned earlier, the tridiagonal structure of the 
Jacobian clearly indicates that while performing the Gaussian 
elimination, only four blocks - the B block being reduced, the 
adjacent C block and the A and B blocks of the next row - are to 
be operated upon at a time in addition to the right hand side 
containing the discrepancy functions. The original A, B and C 
matrices need not be stored, only the four blocks which are being 



operated upon can be generated at a time and the elements of the 
upper-triangular matrix formed after the Gaussian elimination can 
be stored along with the right hand side. 

In the proposed algorithm, the A, B and C blocks and right 
hand side elements which are operated upon at a particular step 
are put in locations g(i,k) where i = l,4c+2 and k = l,4C+3. Then 
the Gaussian elimination from that step is performed. The final 
matrix elements from the j-th row (obtained from & Cj blocks) 
are stored at locations h^(i,k) while the right hand side elements 
are stored at locations d^ (k) . Then the elements of the next row 
of blocks are shifted from their present locations to locations 
where the previous blocks of elements were kept , elements of the 
next generated blocks are put in the last (2C+1 x 4C+3) locations 
of g(i,k). This is continued upto the Nth stage to finally obtain 
the upper-triangular matrix which is then solved by the back 
substitution. 

The elements of blocks B^ and at the jth stage are g(m,k) 
where m goes from 1 to 2C+1 and k goes from 1 to 4C+2. Similarly 
the elements of and kept at locations g(m,k) , m 
going from 2C+2 to 4C+2 and k having the same ranges as before. 
The Fj and Fj_|_ ^vector elements are at locations g(m,4C+3) where m 
= l,4C+2. 

The steps of the Gaussian elimination : 

The forward elimination : 

To perform the forward elimination, a few basic operations 
are defined which can be used with different indices in the 



various steps of the elimination. These operations are as follows 


Operation 1 : 

g(m,k) < g(m,k)/g(in,m) 

Operation 2 : 

g(m,k) < g(iti,k) - g(l,n) 

Operation 3. : 

g(m,k) < g(m,k) + g(l,n) 

Operation A : 

g(m,k) < g(m,k) + g(m,l) g(l,k) 

Operation 5 : 

g(m,k) < g(in,k) - g(m,l) g(l,k) 

Operation 6 : 

h^(m,k) < g(m,n) (For m= 1,2C+1; k= 1,C+1; 

n= C+1,2C+1 

& m= 1,2C+1; k= C+2,2C+2; 
n= 3C+2,4C+2) 


Operation 7 : 

(m) < g(m,k) (For m = 1,2C+1 ; k = 4C+3) 

Operation ^ : 

g(m,k) < g(2C+l+m,k) 

(For m = 1,2C+1 ; k = 1,2C+1 & 4C+3) 


The steps of the Gaussian elimination now become : 


Block (j = 1) ; 

Pivots (1,1) to (C,C) (i = 1 to C) 

Operation 2 — m = C+1 to 2C+1 & 4C+2; k = C+i ; 1 = m ; n = i 
Operation 3 — m = C+1 to 2C+1 ; k = 3C+l+i ; 1 = m ; n = i 



4C+3 


in = 2C+l+i ; k = 4C+3 ; 1 = m ; n = 

m = 4C+2 ; k = 3C+l+i ; 1 = m ; n = i 

Operation 4 — m = C+l to 2C+1 & 4C+2 ; k = 4C+3 ; 1 = i 


Pivots (C+l, C+l) 
Operation 1 — 
Operation 5 — 
Operation 3 — 


to (2C,2C) (i = 1 to C) 

m = C+i ; k = C+l+i to 2C+1 & 3C+2 to 4C+1 
m = C+l+i to 2C+2 & 4C+2 ; k = as above ; 1 
m = 2C+l+i ; k = as above ; 1 = C+i ; n = k 


Pivot (2C+1,2C+1) 

Operation 1 — m = 2C+1 ; k = 3C+2 to 4C+3 

Operation 5 — m = 2C+2 to 3C+1 & 4C+2 ; k = as above ; 1 = 

Operations 6, 7 and 8 are also done. 

Blocks to B,^, ^ (j = 2 to N-1) : 

Pivots (1,1) to (C,C) (i = 1 to C) 

Operation 4 — m = C+l to 2C+1 ; k = C+l to 2C+1 & 4C+3 ; 1 

m = 2 C+l ; k = as above ; 1 = i 

Operation 3 — m = 2C+l+i ; k = as above ; 1 = m ; n = i 

m = C+l to 2C+1 ? k = 3C+l+i ; 1 = m ; n = i 
m = 4C+2 ; k = 3C+l+i ; 1 = m ; n = i 

Pivots (C+l, C+l) to (2C,2C) (i = 1 to C) 

Operation 1 — m = C+i ; k = C+l+i to 2C+1 & 3C+2 to 4C+1 

Operation 5 — m = C+l+i to 3C+1 & 4C+2 ; k = as above ; 1 : 

Operation 3 — m = 2C+l+i ; k = as above ; 1 = C+i ; n = k 

Pivot (2C+1,2C+1) 


& 4C+3 
= C+i 


2C+1 


= i 


& 4C+3 
C+i 



Operation 1 — m = 2C+l ; k = 3C+2 to 4C+3 

Operation 5 — m = 20+2 to 3C+1 & 4C+2 ; k = as above ; 1 = 2C+1 
Operations 6, 7 and 8 are also to be done. 

Block (j = N) : 

Operations are similar to those with blocks B^ to B^^, but the 
steps involving m & k > 2C+1, except k = 40+3, are not to be 

performed. 


The back substitution : 

The basic operations for the back substitution are : 
Operation 1 : 

d^ (k) < d^ (k) - Y. (m) AX(l,i) 

m 


Operation 2 : 

d^ (k) < d^ (k) - AX(j + l,0+k) 

Operation 3 : 

d^ (k) < d^ (k) + AX(j+l,C+k) 

Operation ± : 

d^ (k) < d^ (k) 

Operation 5 : 

f AX(j,k) = d^(k) 

Now the back substitution steps can be written in terms of the 
above operations as follows : 


Last (Nth) B block (j = N) 

AX(j,2C+l) : k = 2C+1 

Operation 6 

AX(j,2C) to AX(j,C+l) : k = 2C to C+1 



Operation 1 — m = C+1 to k+l-C, 1 = N, i = C+m 

Operation 6 

AX(j,C) to AX(j,l) : k = C to 1 

Operation 1 — m = 1 to C+l, 1 = N, i = C+iti 

Operation 4 
Operation 6 

(N-l)th to second B blocks (j = N-1,2) 

AX(j,2C+l) : k = 2C+1 

Operation 1 — m = C+2 to 2C+2 , 1 = j+1, i = m-l 

Operation 6 

AX(j,2C) to AX(j,C+l) : k = 20 to C+l 

Operation 1 — m = C+l to k-C+1, 1 = j, i = C+m 
Operation 1 — m = C+l to 2C, 1 = j+1, i = m 
Operation 6 

AX(j,C) to AX(j,l) : k = C to 1 

Operation 1 — m = 1 to C+l, 1 = j, i = c+m 

Operation 2 
Operation 4 
Operation 6 

First B block (j = 1) 

AX(j,2C+l) : k = 2C+1 

Operation 1 — m = C+l to 2C+1, 1 = j+1, i = m 

Operation 6 

AX(j,2C) to AX(j,C+l) : k = 2C to C+l 

Operations are same as in the case of k = 2C to C+l 



for 2nd to (N-l)th rows 
AX(j,C) to AX(j,l) ; k = C to 1 

Operation 2 
Operation 3 
Operation 4 
Operation 6. 

3.3 Initial Guesses , Criterion of convergence and Mapping 

Functions : 

To ensure convergence with the Newton-Raphson method ,a good 
initial guess of the variables not too far away from the solution 
is required. Different initialisation schemes have been proposed 
in literature. Two of them are mentioned here. 

1) Wayburn-Seader algorithm : In this method, proposed by Wayburn 
and Seader [10], stage temperatures are linearly interpolated 
between two or more guessed values .Usually the first and the last 
temperatures are guessed. The initial guesses of total flow rates 
are computed by considering constant molar overflow. The components 
flow rates are then determined using bubble point calculationa by 
the tridiagonal matrix method. 

2) Initialisation with feed conditions :This is an extremely 
simple method of generating initial guess values. All the stage 
temperatures are set equal to the feed temperature. The total 
liquid and vapor flow rates are calculated by considering constant 
molar flow. The stream compositions are set equal to the feed 
composition. 



The criterion of convergence can be selected from any of the 
various available methods. In the present work , the convergence 
criterion,!^ has been chosen following the recommendation of Henly 
and Seader [3] as : 


2 -10 

= N(2c+1) (1 F.^ ) 10 (3.3.1) 

1 ^ 

Where Fj =total feed to stage j. According to them, this convergence 
criterion will provide results accurate to four or more 
significant digits .The Newton-Raphson iterations are stopped when 
the sum of the squares of the errors are less than or equal to ■il,t 
is , 

j=l 


To avoid the possibilty of flow rates or temperatures becoming 
negative during iteration, a mapping function proposed by 
lenly-Seader [3] has been used. When negative flow rates and 
-emperatures are encounterd, the vai^iables are computed from 
the following relation 


^Itl 


exp 


- -i — 



istead of equation (3.1.12). 



3.4- Outline of Computer Program : 

Computer progarms have been written in FORTRAN 77 to implement the 
modified Gaussian elimination for single columns using the 
Newton-Raphson technique . The programs are divided into subroutines 
to perform the various steps. 

The thermodynamic properties and feed specifications are read from 
separate input files . Parameters like reflux ratio, bottom flowrate 
or fraction of side streams are read from console. 

There are provisions to calculate thermodynamic properties in 
several ways. Like K can be calculated also by :l)As a temperature 
polynomial [ 12 ] 2) From the Wagner equation [6] 3) By gama-phi approa 

apart from the general j Equation of State approach. The progarm 
allows the calculation of enthalpies in two ways :l)As a 
temperature polynomial 2) Using C as a temperature polynomial. 

It 

The partial derivatives , which are the elements of the 
Jacobian , are calculated from analytical expressions in the 
program. The value of the fractional correction in the 
Newton-Raphson method can be calculated from either Broyden's 
formula or can be given by user at each iteration. 

The problems soTved using these progarms are reported in 


chapter 5 . 



Chapter : 4 


HIGH PRESSURE PHASE EQUILIBRIA : EQUATION OF STATE APPROACH TO 
CALCULATE K. 


Distillation is probably the most important separation process 
used in the chemical industry , where the coexisting phases are 
vapor and liquid. Hence reliable vapor liquid equilibrium data are 
required for the design of distillation columns. As it is not 
possible to get experimental VLE data for all multicomponent 
systems, so it is necessary to estimate the VLE data. 

At low to moderate pressures ,the 
Gama-Phi approach can be used to get K.In this approach the 
fugacity coeficient is used to describe the vapor phase and the 
definition of activity coefficient is used to describe the liquid 
phase. At low pressures one can approximate the fugacity 
coefficient as equal to one in vapor phase. 

Vapor liquid equilibria at high pressure are 
conveniently calculated by using an equation of state applicable 
to both the phases. The major advantage of this approach is it's 
applicability over a wide range of pressures . High pressure vapor 
phase equilibria are often expressed in terms of K factors. The K 


factor for any component i is defined as 


Y. 

K. = -i- 

1 X . 

1 


(4.1) 


where y^^is the mole fraction of component i in the vapor phase and 


x^ is it's mole fraction in liquid phase . Equation (4.1) is rewrit 
in terms of fugacity coefficients. The condition for phase 
equilibrium is 



where is fugacity of component in vapor and liquid phases 

respectively . Since for each phase fugacity coefficient is defined 
by 


and 


V 


<P. 


=v 


YiP 


*^1 x^p 


From equation 4.1 it follows that 



4.2 


The (p^ and ipT are to be calculated from some equation of 
state, which in our case is the Peng -Robinson equation of state, 
basic thermodynamic relation to express fugacity coefficient is 
given by 


In (p^= 


-- 

RT -^v 


(£E ) 

^an^^T,V,nj 


RT_ dv -In Z 
V ^ 


4.3 


where the compressibility factor Z is given by 



Mixing Rules : 

The equations of state are generally expressed for pure components 



first and then extended to mixtures. The expressions require the 
called mixing rules. These have the purpose of representing a 
property of a mixture in terms of the composition and properties 
of the pure components. 

Typically a mixing rule expresses a 
mixture parameter in terms of compositions and pure 

components' parameters according to the relation 


‘m 


E I 
i j 




Q. 


13 


4.4 


In this equation the mole fraction yj_yj ®ay apply to a liquid or 
vapor phase /Qjj involves properties of pure components i and 
j. The difficulty arises when interaction terms are 
considered. Rules used to find are called combining rules. 

Van der Wall initially proposed a mixing rule of the type 


"^ij 


Q. .+Q. . 
~ ~ 2 


Q 


ij 


(QiiQj j ) 


1/2 


These are very simple mixing rules. Better mixing rules are 
available now. The one which has been used in the present work is 
Wong and Sandler' mixing rule [11]. 


Outline of Wong and Sandler's mixinf rule with P-R equation of 
state: 

The Peng -Robinson equation of state is given by 



p= 


RT_ 

V-b 


+2bV-b^ 


where 


0.07780 RT 
c 

c 


.45724 R^T^ 
c 


and a= [l+fw(l-T^/^] ^ 

P ^ 

c 

where fw = .37464 +1.54226cj - .26992w^ 


Based on the statistical thermodynamics the following mixing 

rule is given by Wong and Sandler [11]. 

The expressions for b and a are 

mm 



Q 

(1-D) 


^m (i-D)^ 


with Q and D defined as 


and 


Q=I I x.x. (b-a/RT).. 

i j j J-j 


^ ? ^i b^RT' 



— 00 


CRT 


The fugacity coefficient is calculated from equation 4.3 

For the P-R equation of state with any arbitrary set of mixing 

rules for a^ and b^ we have (reference can be made Wong and 



Sandler [11]) . 
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The partial derivatives of a and b are : 
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with the partial derivatives od Q and D given by 
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hough any free energy model can be used , here UNIQUAC has been 
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:n this case ,the partial derivatives of A ^ /RT with respect to 

:he mole number of each species , which is the logarithm of the 
specific activity coefficient , is given by 


cp . 1? . ^ . N 

Lnri- * I 


N 


q.ln(E i^.r..) + 

j 31^ 


- 


N 1?.T. . 

g r 

^iV N 

^ E 1^1 ■'^1 • 

^ k k: 


where 


<P: 


riXi 


N 

E r^x 

k 


and 


k^k 


T? .=- 
1 


qi^i 


N 




here r and q are pure component parameters and the coordinatio 
number z=10. 
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Chapter : 5 
Results and Discussion 


The modified Gaussian elimination for single column was 
implemented on HP 9000 Supermini computer. The code was written in 
FORTRAN 77. The single column program has been tested with 
distillation problems. The details are presented below. 

TEST OF THE PROGARM : 

PROBLEM 1 (Naphtali and Sandholm 1971 ) 

This problem has been taken from Naphtali and Sandholm paper 
[5]. The components to be separated were n-hexane and 
ethanol. The data for enthalpy calculation, critical properties were 
taken from Prausnitz and Reid [6] 

The UNIQUACparameters were taken from Prausnitz et.al. [6] .The 

feed was vapor at 200°F , consisting of 50 moles of each 

component . The column contained ten plates with the feed entering 

the sixth. The reflux ratio was specified as 5.0 and the column 

was to recover 70% of the feed in the bottom product. This problem 

converged in seven iterations toan error criterion of the order of 
—5 

10 .The results obtained were almost exactly matching with that 
of Naphtali and Sandholm' s. 

The nature of output file is shown on the next page. 
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• i mu 1 411' 1 0 n of -i d t ii i i 1 i * t i o o c. o 1 u m n 


number 

of 

stages 

X 

1 0 

number 

of 

components 

- 

£ 

number 

of 

reboilers 

X 

1 

type of 

T' e b 0 i 1 e r 


par t i ai 

total bottoms flow rate 

X 
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of 

condensers 

= 

1 

type of 
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reflux 
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a 
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t 

£ 

} 
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C- 
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:$ 
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q 
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Problem * 1 - • this problem, a distillation column with a large 

number of stages has been simulated. A 72 stage column separates a 
binary mixture into nearly pure products [3]. 

The feed is a saturated liquid at the feed pressure of 9.4 atm. 
and temperature of 363K. 

The K values are calculated from the Wagner equation. The 
vapour enthalpies are calculated in terms of , taken as a 

polynomial function of temperature. The latent heats are 
subtracted from the vapour enthalpies for computation of liquid 
enthalpies. The details of the output file are presented in the 
next page. 

The problem converged in 11 iterations using the 
initialisation procedure with feed conditions as described in 
section 3.5. 0 was taken to be 1.0 for all the steps in this 

case. 

Problem : 3 

This problem deals with the simulation of a 
distillation column .A superheated vapour mixture consisting of 
five components is separated by a distillation column having an 
overhead partial condenser and a partial reboiler . The column 
consisting of thirteen stages and is operated at 400 psia.Feed 
enters the column at seventh stage at 105 °F .The feed consists of 
a mixture of to compounds . 

The initial stage temperatures 
are set equal to the feed temperature . The total liquid and 
vapour flow rates are calculated by assuming equimolal flow .The 
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6 f ’ 

. 1 6554E-02 

1 25 . 65 

63 

. 1 6566E -02 

125.65 

5 4 

. 1 650 1 E-02 

125 . 65 

65 

. 1 65 6 02-0 2 

1 ES . 65 

6 6 

. 1 65665-02 

125.65 

5 7 

. 1 C564E~02 

125.65 

56 

. 1 6559E-02 

1 25 . 65 

6 9 

. 1 05212-02 

1 25 . 6S 

TO 

. 1 65092-02 

125.65 

71 

. 1 6504E-08 

1 25 . 85 

72 

. 1 66 06E-02 

1 £5 . 85 


liquid phase component flow rates 


stage nos. component no. 


1 

306 . 29 

37 . 324 

£ 

1 07 4 1 

161 .58 

3 

27. SI 4 

231.19 

4 

14.SS9 

246 . 96 

S 

12. 707 

249 .40 

6 

1 2 . 44? 

249 . 75 

7 

12.410 

<?4 9 . 3 0 

8 

12 405 

2 4 9 81 

9 

12 . 405 

249 , 81 

1 0 

12,405 

£49 .81 

1 1 

12.405 

249.81 

12 

12,405 

249.81 

13 

12.405 

249.81 

1 4 

1 2 .405 

£49 .81 

1 5 

12.405 

£49.81 

1 6 

12.405 

249 .81 

1 7 

12.405 

249.81 

18 

1 2 .405 

£49.81 



t 9 

\ a. , 4 OS 

2 4 9 8 1 

c 0 

1 £• , 4 05 

2 4 9 , a 1 

?.) 

12 ’105 

249 . 81 

P j> 

1 2 . 4 05 

249 . 81 

£‘3 

1 2 . 4 0 5 

249 . 81 


1 £ . 4 05 

249.81 

P 

1 2 . 405 

249 . 81 

H6 

1 2 . 405 

249 - 81 

£7 

1 2 . 405 

249 .8? 


1 2 . 4 OS 

249.81 

£9 

12.405 

249 .81 

30 

12.405 

249 . 81 

3 1 

1 c* .405 

249.81 

32 

12.40 5 

249 . 81 

33 

1 c: . 4 0 5 

£49 . 81 

3^ 

12.405 

249 . 31 

35 

12.405 

249 .81 

36 

12.405 

£49 .81 

3T 

12.405 

249 . 81 

38 

1 2 40S 

249 . 81 

39 

12 4 05 

249 . 81 

4 0 

1 2 , 4 0 5 

24 9 , 8 1 

4 1 

1 1 .402 

229 . 6 1 

4 a 

£ . 7 032 

242 . 5 1 

4 3 

f.- O JO p 

245 . 88 

44 

. 1 2T4 9 

246 . 63 

43 

. 2T£96e-01 

£46 , 80 

4 6 

.717912-02 

£46 . 83 

^7 

. 394882-02 

£46 . 84 

48 

337302-02 

£46.84 

4 9 

. 325765-02 

246 . 84 

50 

. 32321 2-02 

246 , 84 

St 

. 318022-02 

246 . 64 

52 

. 31 585E-02 

246 . 54 

53 

. 31823E-0E 

£46 . 84 

54 

. 31805E-02 

£46 . 84 

55 

. 31 790E-02 

246 . 84 

56 

.31801E-02 

246 . 84 

5T 

. .3t810E-02 

246 . 84 

58 

. 31 74SB-02 

246 . 84 

59 

. 31 726E-02 

246 . 84 

60 

. 31 724E-02 

246 . 84 

6 1 

. 31 T29E-02 

246 . 84 

62 

. 31 735E-02 

246 . 84 

63 

. 31 744E-0E 

246 . 84 

64 

. 31 T48E-02 

246 . 84 

65 

. 31 T49E-02 

246 .84 

66 

. 31 74SE-02 

246 . 84 

67 

. 317405-02 

246 . 84 

58 

. 31 732E-02 

246 . 84 

6 9 

. 31 72TE-0E 

246 . 84 

70 

. 31T24E-02 

246 . 84 

71 

.317455-02 

246 .84 

72 

.1S9.75E-02 

121.00 


final vapour, 1 iquid rate and temperature values 


t 

i 15.200 

345.6153 

266.7 

2 

460 .815 

£68 , 9984 

288'. 0 

3 

384 . 1 98 

£58.7038 

317.8 

4 

373.904 

£61 ,5161 

326.7 

5 

376.715 

£62.1107 

328 . 1 

6 

377,31 1 

£62 . 1 977 

328.3 

T 

3?T .398 

tC 6 2 . 2 0 9 8 

328 . 3 



6 

3 T 7 . 4 1 0 

86 8 

.8115 

3£8 . 3 

9 

377.412 

868 

2117 

388.3 

1 0 

377. 41£ 

862 

21 18 

388.3 

t 1 

377.418 

£6 2 

2118 

388. 3 

1 £ 

377,418 

868 

8118 

388 . 3 

1 3 

377.418 

£6£ 

8118 

328. 3 

1 ^ 

377.418 

86 8 

8118 

388. 3 

1 B 

377 .418 

868 . 

8118 

328.3 

1 6 

3TT . 41 £ 

£68 . 

2118 

328.3 

1 T 

377.418 

8G8 , 

£113 

388. 3 

1 8 

377 .418 

868 . 

8118 

328.3 

1 

377.418 

868 . 

81 1 B 

388 . 3 

£0 

377.418 

868 . 

8113 

328 . 3 

£1 

377.418 

868 . 

2 1 1 S 

3£8 . 3 

££ 

3T7 .418 

£68 

8118 

388 . 3 

£3 

377.418 

868 

8 1 1 8 

328 . 3 

£4 

377.418 

868 

£118 

328 . 3 

£5 

377.418 

86 8 

2118 

328 . 3 

£ 6 

377 418 

£68 

£118 

328 . 3 

£T 

3 T 7 4 1 8 

8 6 c 

8 116 

328 . 3 

£8 

377418 

£62 

8 116 

388 . 3 

£9 

3 7 7 4 18 

SCsa 

C-’ 1 1 8 

388 . 3 

30 

3 r T 4 1 8 

8 68 

.2118 

388 . 3 

31 

377 418 

8 6 8 

£ 1 1 8 

328 . 3 

32 

377.418 

£ 6 8 

£ 11 3 

328 . 3 

3 3 

377418 

268 

£118 

328 . 3 

34 

377.418 

868 

£113 

328.3 

.35 

377 .418 

£68 

£ 1 1 8 

323 . 3 

36 

377. 418 

£ 6 £ 

2118 

328.3 

3 7 

377.418 

868 

2118 

328.3 

38 

377 . 418 

£68 

8118 

328 . 3 

39 

377.418 

868 

8118 

3^8.3 

40 

377 .418 

£68 

£118 

328.3 

41 

377.418 

£4 1 

0 110 

328.3 

4£ 

120.018 

845 

8148 

336 . 6 

4 3 

184.881 

846 

•4 738 

338 . T 

44 

1 £5 . 460 

84 6 . 

7680 

339.2 

45 

125.768 

£46 . 

8 84 0 

339.3 

46 

185.830 

£46 . 

8385 

339 .3 

47 

1 85 . 843 

£4 6 

8486 

339 .3 

48 

185. 847 

£46 . 

8437 

339.3 

49 

1 £5 . 848 

£4 6 . 

8438 

339 .3 

50 

185. 648 

£4 6 

8437 

339.3 

SI 

185.848 

£4 6 . 

84 34 

339.3 

sa 

185.847 

246 . 

8 4 35 

339.3 

53 

125.847 

246 . 

8427 

339.3 

54 

125 . 847 

84 6 

8487 

339.3 

55 

185.847 

24 6 . 

S434 

339.3 

56 

185 . 847 

24 6 . 

8433 

339.3 

57 

185.847 

246 

8439 

339.3 

58 

185.848 

£46 

8434 

339.3 

ro 

125 847 

P 4 6 

r 4 3*^ 

339 . 1 

60 

I .',2 . i\4 f 

. H. 

fi .} i 

1 .1 0 . S 

61 

185 . 848 

£ 4 6 

8443 

339.3 

6£ 

185.848 

£46 

8441 

339.3 

63 

185 848 

246 

8447 

339.3 

64 

185.849 

£4 6 

8456 

339.3 

65 

1 £5 . 850 

£46 

8456 

339.3 

66 

185.850 

246 

8457 ' 

339.3 

67 

185.850 

246 

8452 

339.3 

66 

185.849 

24 6 

8 456 

339.3 

69 

185.849 

246 

84 65 

339.3 

7 0 

185.850 

246 

8452 

339.3 

71 

12S.849 

846 

8451 

339.3 

78 

1 85 649 

181 

0 016 

339.3 


Temperature 1 


stage 


V 


y 



1 

£66.66 

2 

£9 9 0 3 

3 

3 1? . e 0 

4 

326 . 69 

5 

, 1 0 

6 

3 c 8 . 3 0 

7 

32:0 , 3 3 

8 

328 . 33 

9 

3.?!? 3 3 

1 0 

3 £ 9 3 y 

1 1 

3 «;■' 9 3 3 

1 £ 

3c 8 3 3 

1 3 

32S . 33 

1 4 

3£8 , 33 

1 5 

3r-?. , 3 3 

1 6 

328 . 33 

1 7 

328 . 3 3 

1 8 

3 i;! 8 ,3 3 

1 9 

328 7 3 

2 0 

3c 9 3 3 

£1 

32 8 3 3 

22 

328 . 33 

£3 

3c 8 . 33 

2 4 

3 £ 8 3 3 

2 5 

328 , 33 

26 

328 , 33 

2 7 

328 . 33 

£8 

328 33 

£9 

388 . 33 

30 

328 33 

31 

328 33 

3£ 

328. .3.3 

3 3 

328 . 33 

34 

388 3 3 

35 

3£8. .3 3 

36 

328 . 33 

37 

3es . .33 

38 

3£8 . 33 

39 

3r?8 , 3 3 

40 

328 . 3.3 

4 1 

388 . 33 

42 

336.58 ' 

43 

338 . 7£ 

44 

339 . £0 

45 

339 . 30 

46 

339 . 33 

47 

339 . 33 

48 

339 . 33 

49 

339 . 33 

so 

339 . 33 

51 

339 . 33 

52 

339 . 33 

S3 

339 . 33 

54 

339 . 33 

55 

339.33 

56 

339 . 33 

57 

339 . 33 

58 

339.33 

59 

339 . 33 

6 0 

339 . 33 

61 

339 . 33 

6£ 

339 . 33 

63 

339.33 

64 

339 . 33 

65 

339 , 33 

66 

339 , 33 


3-15 . G£ 
£69 00 

ess ?o 

£ 61 .SB 
eea. i t 
£62. £0 
£6c' . ei 

£6 e . £ 1 
£6£ . £1 
£6£. £1 
£68 £1 
£68 . £1 
e6£ . £1 
£6£ . £1 
£62 . £1 
eee . £i 
e6£ £1 

e G£ . £ t 
eee , £i 
eee . cm 
e 6 e . £ t 
e 6 £ . £ I 

e 6 e £ 1 
£60 £1 
c ! 6 e , £ 1 
H 6 £ . £ 1 
£ 68 . £1 
£ 68 , £ 1 
e 6 £.£l 
2 6 £ . £ 1 
262 . £1 
£ 62 .21 
£ 62 .21 
862 21 
262 .21 
262 21 
262.21 
262 .21 
262.21 
2 6 2 . 2 1 
241 .01 
245 .21 
246 . 47 
246 . 76 
246 . 82 
246 . 84 
246 . 84 
246 . 84 
246 . 84 
246 .84 
246 . 84 
246 . 84 
246 . 84 

246 . 84 

246 . 84 

246 . 84 
246 . 84 
£46 . 84 
246 .84 
246 . 84 
246 . 84 
246 . 84 

246 . 84 

246 . 85 
£46 .85 
246.85 


M 5.20 
460.82 
384 .20 
373 . 90 
376 . 72 
377,31 

377 .40 

377 .41 
377 . 41 
377. 41 
377 . 41 
377.41 
377. 41 
377 . 41 
377. 41 
377. 41 
377. 41 
377 41 
377. 41 
3 77 4 1 
377.41 
377. 41 
377 . 41 
377.41 
377 . 41 
377 4 1 
377 41 
377 .41 
377 .41 
377 .41 
377 .41 
377.41 
377 . 4 1 
377 41 
377.41 
377.41 
377 . 4 } 
377. 41 
377.41 
377 . 41 
377.41 
120.02 
124.22 
125.48 
125.77 

125.83 

125.84 

125.85 
125.85 
125.85 
125.85 
125.85 
1 25 . 85 
125.85 
125.85 
125.85 
125.85 
125.85 
125.85 
1 25 . 85 
125. 85 
125.85 

1 25 . 85 
125.85 
125.85 
125.85 


S9 eE'f 0 0 
399E-f 00 
t 0 65+0 0 

- ^855-01 
. <1755-01 
■ 473E-01 
- <173E-“01 
. 4735-01 
. 473E-01 
. 473E-01 
. 4 7 3E- 0 1 
473E-01 
> 473E-0t 
■ 47.3E-01 
. 473E-01 
. 473E-01 

. 4 7 35-01 
. <*7 35- 01 
4 7 3E--01 
. 4T3E- 01 
. 4T3E-01 
. 4 7 3E- 0 1 
.4732-01 
. 473E-01 
• 473E-01 
4735-01 
. 473E.-01 
473E-01 
. 4 7 3E-0 1 
. 473E~01 
. 4 7 35” 0 7 
. ^73E~ 0 1 
^?3E • 01 
4 r3E-01 
. 4 73E-01 
. 473E-01 
. 473E-01 
. 4 T3E.. 0 1 
. 4 73E- 0 1 
. 1 1 OE-01 
. 241 e:-02 

■ SI 7E-03 
. 1 1 IE-03 
.c!91E-04 
. 1 60E-04 

■ I 37E-04 
. 1 38E-04 
. 1 31 E-04 
. 1 29E-04 
. 1 2SE-04 

. 1 2 9E-0 4 
. 1 29E-04 
. 1 29E-04 
. 1 c:9E" 0 4 
. 1 04 
. I E9E-04 
. 1 29E-04 
. 129E-04 
. I29E-04 
. I e9E- 04 
. 129E-04 
- 1 29E-04 
. 1 29E- 04 
. 1 29E-04 


. 99SE+00 
- 918E+00 
.578E+00 
. 380E+00 
.343E+00 
.33TE+00 
. 337E+00 
.337E+00 
.33TE+00 
. 337E+00 
.337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+ao 
. 337E+00 
. 337E+0a 
.337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
.337E+00 
. 337E+00 
- 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
- 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 337E+00 
. 950E-01 
.218E-01 
.472E~02 
• 1 01E-02 
. 217E-03 
. 487E-04 
. 1 99E-04 
. 152E-04 
. t43E-04 
. 1 40E-’04 
. 122E-04 
. 1 33E~04 
. 1 34E~04 
1 33E-04 
133E-04 
1 34E-04 
134E-04 
t32E-04 
1 31 E-04 
I 31E-04 
1 32E-04 
1 32E-04 
1 32e-04 
1 32E-04 
1 32E-04 



ST 

130 . 33 

£46 . 85 

1 £5 . 85 

. 1 295-04 

. 1 32E-04 

SB 

139 , 33 

£46 . as 

T£5. 65 

, f a'JE-94 

. 1 3£E“04 

S3 

339 . 33 

£46 , BB 

leS . 85 

. 1 c!9E'-04 

. 1 31E-04 

TO 

339 . 33 

£46 , 85 

125.85 

t 29E-04 

. 13t E“04 

71 

3 3 9.33 

£46 . 85 

125.65 

, t 2 9E- 0 4 

. 131 E-04 

72 

339 , 33 

121 .00 

1 25 . 85 

1 32E“-04 

. 13£E-04 


condenf^er 

heat duty == 

.S08654E + 0 ? 




r e b o 1 1 e r 

heat duty » 

. 3697435+07 




number of 

iterations = 

1 1 





compositions of each liquid and vapour stream are set equal to 
the composition of the feed . 



fnoBLGtM-S 

simulation of a d I f t i 1 1 at i on column 
n u m b r o f ?•* t a g e *5 •- 13 

number of c o n.> p o n e n t h-: 5 

number of re boilers = 1 

type of reboiler ; partial 

total bottoms flow race = £80,00000 

number of condensers = 1 

type of condenser : partial 

ref 1 u X rat io ~ 1,890 


results of s i mu 1 a t i o n 


error criterion - .91S300E-08 

vapour phase component flow rates 
stage nos, component no. 




' 

£ 

3 

4 

5 

1 


1 60 . 00 

358 . 92 

1 0815 

- T6531 E-02 

. 1 9824E-02 

£ 


£25 . 1 £ 

1 265 . 3 

1 2 . 3 T 9 

. 290435-01 

. 66728E-02 

3 


189.12 

1278.5 

35 . 376 

. 3491 3E-01 

-70T68E-02 

4 


1 82 . 95 

1200 . 4 

87.018 

. 6T525E-01 

.731E1E-02 

5 


180.60 

1 057 . 5 

1 79 . TT 

. 36424 

-9065tE-02 

6 


178.7^ 

887 . 99 

297 . 88 

2.0406 

. 47394E-01 

7 


ITT. 25 

742 . 95 

397 . 54 

8 . 7394 

.57027 

8 


1 5 . 674 

265 .37 

190 73 

4 . 0287 

. 26479 

9 


3 . 721 4 

236 . 1 8 

246 20 

4.960 1 

. 32603 

1 0 


. 83093 

186.57 

302 47 

5 . 9003 

.38451 

1 1 


. 20324 

133.75 

358 65 

T. 1504 

.44980 

1 2 


, 85631 E-01 

86 . 285 

407 . 83 

9 . 7464 

.58438 

I 3 


. 559372-01 

47 , 974 

437 IT 

1 7 . 686 

1 .4334 

liquid 

phase 

component 

1 

flow rates 




Stage 

nos . 

component no. 

_ 





1 

2 

3 

4 

5 

1 


65.121 

906.38 

1 t . 29T 

.21T20E-01 

.47174E-02 

2 


29.124 

919.58 

34 , 294 

279585-01 

. 51 408 E- 02 

3 


22 . 953 

841 44 

85 . 936 

. 61 T65E-01 

.53T72E-02 

4 


20 .600 

698.65 

1 78 6 9 

. 36388 

.71595E-02 

5 


1 8 . 739 

529 . OT 

£96 . 80 

2 . 0406 

. 47060E-01 

6 


1 T . £55 

384 . 03 

396 . 46 

8 . 7394 

.57027 

T 


1 5 . 674 

2T6 .45 

429 . 64 

29 . 029 

5.2643 

8 


3 . T21 8 

247 . 26 

485 . 1 £ 

. 29.960 

5 . 3260 

9 


. 63456 

197.65 

541 39 

30 . 900 

5 . 3845 

1 0 


. 22574 

144.83 

597 57 

32 150 

5 . 4498 



t 1 

. * ! 'SB'S 

9 7 . 36 6 

6 *\ 6 ? S 

34 . 746 

5 . 5844 

1 a 

. e3988E-“0t 

59 . OSS 

676.09 

4£ . 686 

6.4334 

1 3 

. E87aiE-01 

11.081 

£38 . 92 

25,000 

5. 0000 


final 

vapour , 1 i qui d rate 

and temperature values 




1 

520 01 0 


982 . 8267 

11.66 


e 

1502 . 54 


983 . 0302 

28 .40 


3 

1503 . 04 


950 .3916 

34, TT 


4 

1 470 . 40 


898 .3143 

43.80 


5 

1416.32 


846 . 6976 

57.77 


6 

1366.70 


807 . 0547 

73,65 


7 

1 327 . 05 


75r> . 0645 

88 . 23 


8 

476 . 064 


77 T . 3887 

105.7 


9 

491 .339 


776 .1581 

1 IT. 7 


1 0 

496.155 


780 . 2266 

128.7 


1 1 

500 . 205 


784 . 5571 

138.6 


1 £ 

504 . 529 


784 . 3445 

147.5 


1 3 

504 .317 


£80 . 0283 

157.4 

St age 

Temperat ure 

1 

V 

X 

y 

I 

n . 6 6 

982,63 

520 . 01 

. ersiE-oi 

.30854-00 

E* 

28.40 

983 . 03 

1502. 84 

. 2962-0 1 

. 1 50E + 00 

3 

34 . T7 

950.39 

1503 . 04 

. 242E-0 1 

. 1 26E+00 

4 

43 . 80 

898 . 31 

1 470 . 40 

. 229e--01 

. 1 24E+00 

5 

5T . 7 7 

846.70 

1418.32 

. 2215-01 

. 1 27E+00 

6 

73 65 

507 . 05 

1 366 . TO 

.2145-01 

. 1 31E+00 

7 

68 23 

756 . 06 

1 327 . 05 

. 2075-01 

. 1 34E + 00 

8 

1 05 72 

771 .39 

476 . 06 

, 48^E-02 

.329E-01 

9 

117.75 

776.1 6 

491.39 

. 1085-02 

. 757E-02 

1 0 

128. . 6 6 

780 . 23 

496 . 15 

. 2895-03 

.167E-02 

1 1 

1 38 . 62 

784 . 56 

500.20 

. 1 455-03 

.406E~03 

1 2 

147.48 

784 . 34 

504 . 53 

. 1 075-03 

. 1 70E-03 

1 3 

157.36 

280 . 03 

504. 32 

. 1 035-03 

. 1 TIE-03 


condii»nser heat duty * .474591E+07 
reboiler heat duty ~ .£9£744E+07 
number of iterations = 7 



Conclusions 


The Gauss elimination has been suitably modified to solve 
the Jacobian in case of single columns. The inner sparsity of the 
blocks of the block-tridiagonal Jacobian has been exploited. The 
equation of state has been sucessfully incorporated to take care 
of highly non-ideal systems. The beauty of this approach tc 
calculate equilibrium constants [K] is this that it does not requ 
experimental data. It needs only critical parameters. 

Computer codes have been written to implement the modified 
Gauss elimination technique and equation of state approach. The cod- 
worked well for single distillation columns. 



Nomenclature 


(2C+1) X (2C+1) matrices 

nonzero element at location (i,l) of M 

vector of order m 

number of components 

specific heat 

scalar 

array of order Nx(2C+l) 

enthalpy balance discrepancy function for the 
j-th stage 

. T 

(n+1) order unit vector, [ 0,0,.... 0,1] 

vector of discrepancy functions of order N(2C+1) 

vector of discrepancy functions of order (2C+i) 

flow rate of component i in the feed to j-th stage 

array of order (4C+2) x(4C+3) 

enthalpy of vapour of component i leaving the 

j-th stage 

scalar 

enthalpy of liquid of component i leaving the 
j-th stage 

array to store upper-triangular matrix, order 
Nx(2C+l) x(3C+2) 

array of dimension (2C+1) x(2C+l) to store lower 
off-tridiagonal blocks 

array of dimension (2C+1) x(2C+l) to store upper 
off-tridiagonal blocks 
identity matrix of order CxC 
Jacobian of order N(2C+1) xN (2C+1) 



''j.l 





M 

m, m’ 
N 


n 


0 

0 

p 


c 

cc 



R 


a,b 





i 


equilibrium constant of component i at the 
j-th stage 

component liquid flow rate for component i leaving 
j-th stage 

total flow rate of liquid leaving the j-th stage 
material balance discrepancy function for 
component i leaving stage j 
square matrix of order mxm 
matrices of order NxN 
number of stages 
number of equations 
null column vector of order C 
null matrix of order CxC 
vector of order n 

discrepancy functions for phase equilibrium 
relations for component i leaving stage j 
heat supplied/withdrawn from the j-th stage 
fraction of liquid leaving stage b for stage a 
fraction of liquid stream from (j-l)-th 

stage reaching j-th stage 

fraction of vapour stream from (j+l)-th stage 
reaching j-th stage 

fraction of the vapour side-stream from j-th stage 
fraction of the liquid side-stream from j-th stage 
temperature of the j-th stage 

total flow rate of vapour leaving the j-th stage 
component vapour flow rate for component i leaving 
j-th stage 



X 


a nonzero element 



D/i 


1,1 




V 


j , i 


e 


a 

A 

b 

B 

F 

V 

<P 

r 

r,q,x 


vector of variables of order Nx(2C+l) 
column vector of nonzero elements of order C 
matrix of nonzero elements of order CxC 
vector of variables of order (2C+i) 

desired solution of equation system 
initial guess vector 

mole fraction of component i in the vapour leaving 
the j-th stage 

mole fraction of component i in the vapour leaving 
the j-th stage 
convergence criterion 

Murphree efficiency of component i at j-th stage 
fraction of the Newton-Raphson correction 
equation of state energy parameter 
molar helmholtz free energy 
equation of state excluded papameter 
second virial coefficient 
arbitrary function 
molar volume 
fugacity coefficient 
activity coefficient 
uniquac papameters 
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APPENDIX 


The elements of the A, B and C blocks : 


A. 
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dF. . 

1,1 


dX 


j-l,k 


5M. . 

1 , 1 


= r 


ai 


j-l,k 


T 5 . , 

Lj_i i,k 


5 . , 
i,k 

5 . , 
i,k 


1 , if i = k 
0 , if i k 


5 M. 


1,1 


= 0 


1-1, k 




= 0 


i 3Q, . ag. . 

= 0 ; = 0 ; 


51j_i^i 


dv . . . 

1 - 1,1 


"^j -1 


5E. 


^^j-l,i 


dE. 

= L._^ ^- 1,1 ' — 


= 0 


j-i,i 


aE 






ah 


j-l,i 



aM. . 

D ,1 

= -(l+s.)5. 

V ' 

aM. 

3 

av. - 
D ,k 

J J 

, K. 

81. 
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V K. . 

3 3 t ^ 


_ 7 

"ij,k 




7 ]. .1. . 5K. . 

j , X j , 1 3,1 


5 . , 

l,k 


■q 1 dK 

j,i j,i j,i 

L dv~r 




7]. .1. . 5K. . 

'3,13,1 3,1 


= -d+s )h. . 


-3 = -(l+s.)H. . 

SVj^i 3 ' 


5h. . 

= -( 1 +s^) 11 ^, 

^ i=l ^ ^ 3T . 

3 


an. . 

- ( 1 +s.) z V . . 

^ i=i ^ ' aT j 


c . = 


ar. . 

D /I 


®^j+l,k 


j , i 


aM. . 

3 , 1 


j+i, k 


3 + 1 


/ i 


0 



0 


oQ. . 

: ,1 


ai 


j+l,k 


dQ . . 

dv . , . T 

3+1, k 


( 1 - 7 ]. .) 

^ 3,1^ 




3+1 


3+1 


(''V. ^j+1 + ^ ^V. 

3+1 -' j+1 




dQ. . 
3,1 


aT 


= 0 


3+1 


aE. 

3 


ai 


j+i,k 


^j,3+l^j+l,k ° 


aE. 


av 


j+i,k 




aE. 
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aT 
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•■ Tv. 
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